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Abstract 

Ewald summation is widely used to calculate electrostatic interactions in computer simulations of condensed- 
matter systems. We present an analysis of the errors arising from truncating the infinite real- and Fourier-space 
lattice sums in the Ewald formulation. We derive an optimal choice for the Fourier-space cutoff given a screening 
parameter rj. We find that the number of vectors in Fourier space required to achieve a given accuracy scales with 
?7 3 . The proposed method can be used to determine computationally efficient parameters for Ewald sums, to assess 
the quality of Ewald-sum implementations, and to compare different implementations. 



1 Introduction 

The calculation of electrostatic interactions in computer-simulation studies of condensed-matter systems poses serious 
problems regarding accuracy and efficiency H . These are mainly caused by the infinite range of Coulomb interactions 
in conjunction with the finite size of the samples studied. To avoid system-size effects, periodic boundary conditions 
are employed. The natural description of the electrostatics in this periodic space is obtained by summation of the 
charge interactions over periodically replicated simulation cells. This yields the well-known formula of Ewald Q for 
the Coulomb energy of charges in a lattice. 

Ewald summation is widely used in computer simulations of charged systems and is generally believed to give 
the most accurate description of the electrostatics (with respect to system-size dependence). The Ewald formula 
splits the Coulomb energy into two rapidly converging real- and Fourier-space lattice sums, where their relative 
contributions can be controlled by a parameter r\. However, in numerical implementations one has to apply truncations 
of the two (infinite) lattice sums. In this work, we develop a quantitative description of the errors arising from this 
truncation. We then discuss the choice of cutoff distances in real and Fourier space with respect to numerical accuracy, 
resulting in restrictions regarding computational efficiency of Ewald-sum implementations. In particular, we analyze 
the connection between the real-space screening parameter r\ and the number of Fourier-space vectors required for a 
given accuracy. We also derive an approximate upper limit for an efficient choice of the Fourier-space cutoff. 



2 An accuracy measure for truncated Ewald sums 



The Ewald-summation formula for the Coulomb energy U of m charges q a (with net charge 0) at positions r a can be 
expressed as a sum over pair interactions and self terms, 
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where r a j3 = Yp — r a + n, with the lattice vector n chosen such that r a p is a vector in the unit cell. This result holds 
for a background dielectric constant e = 00 (or vanishing dipole moment of the cell), as discussed in ref. Q. The 
effective pair interaction has the following form, 
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where V is the volume of the box, erfc is the complementary error function, and k = |k|. The two lattice sums 
extend over real and reciprocal (Fourier) space lattice vectors n and k, respectively. In most practical applications, 
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the screening parameter r\ is chosen such that only n = and a few hundred k vectors need to be considered. The 
convergence parameter r\ can be chosen arbitrarily (0 < r\ < oc), as 



(3) 



The choice of r\ gives different weights to the two sums. However, the requirement of numerical accuracy imposes 
some restrictions on the choice of r\. Any truncation of the two lattice sums in Eq. (|J) results in deviations from 
the identity Eq. (|^). This provides us with a measure for the accuracy of a given implementation characterized by 
a screening parameter r\ and two cutoff distances N and K for the real- and Fourier-space lattice sums (N > |n|, 
^>|k|). 

The errors AU of the total energy and Au a of the single-particle energies are weighted sums of numerical errors 
Aip in the electrostatic interaction ^(r) = <p(r) — 1/r, 
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In most practical cases, we expect a considerably smaller error than indicated by the upper bounds of Eqs. (||) and 

t. A detailed analysis of the errors in disordered charge configurations has been presented by Kolafa and Perram 
. In this work, we focus on a model- and configuration-independent measure of the numerical accuracy of truncated 
Ewald sums. This general error analysis provides insight into the influence of the Ewald-sum parameters 77, N, and 
K (and, possibly, a real-space cutoff R c ). We do not analyze effects of partial error cancellation owing to, e.g., charge 
ordering. 

We will study a system of two particles with charges ±1, at positions r = and r. We will restrict our analysis to 
the most widely used cubic cell. Calculations will be done in reduced coordinates with a lattice constant of 1, resulting 
in n assuming integer values and k = 27rn. We define A(rj, N, K) as the maximum deviation from the identity Eq. (||) 
with respect to r vectors in the cell V. A(rj, N, K) will serve as our measure for the numerical accuracy of a truncated 
Ewald sum. From Eq. (|) we obtain 



A(ry, N, K) = max 
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Using the identity Eq. (||) for the full Ewald sum (N, K — » 00), we can invert the sign and sum over the complementary 
n- and k-space regions n > N and k > K, respectively. We approximate the fc-space contributions neglected in Eq. (|^) 
by an integral, 
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noting that A& > 0. This results in an approximate expression for the neglected fc-space contributions, 

A fc (77,X) « *(K/ V ). (8) 
We obtain cr{x) by integration of Eq. (0) replacing 1 — sm{x)/x by its maximum value 7 (~ 1-22), 

a(x) = 27X7r _1 exp(-a; 2 /4) + 277r _:L/2 erfc(x/2) . (9) 
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Typically, 77 is large (77 > 4) and only one simulation cell is considered in r space (N — 0). Regarding the sign of the 
real-space contributions for N = 0, 



A r (r),N = 0,K,r) = 2tt' 1/2 V [exp(-? ? 2 |r + n| 2 ) - exp(-7? 2 7i 2 )] , (10) 



n 

n>0 



A r is found to be positive in extensive numerical tests. Expressed in terms of theta functions using Jacobi's imaginary 
transformation [gj, we conjecture 

> ^V 3 [M™,Q)Mw,q)M™,Q)-*l{o,Q)] > explVte 2 + y 2 + * 2 )] - 1 , (n) 

where $3 (u, q) = l+2X^i<7™ cos(2nit), q = exp(— n 2 /^ 2 ), and (x, y, z) e [— 1/2, 1/2] 3 . The left part of the inequality 
is trivial. However, we could not find a formal proof for the lower bound. The neglected r-space contributions in A 
come mainly from the neighboring image and reach a maximum at the center of the faces of the cube [r = (1/2, 0, 0)]. 
From Eq. (|To|), the real-space contributions to A are estimated as 2-7T -1 / 2 exp(— 7? 2 /4), resulting in an approximate 
expression for A, 

A(r?, N = 0, K ) w a{K/rf) + 2tt" 1/2 exp(-r; 2 /4) . (12) 



For the practically less interesting case of N = and 77 small (77 < 1), Eq. ( |ll| ) can be used to derive an upper bound 
A r < 2-7r~ 1 / 2 [l — exp(— 3?7 2 /4)]. Based on an analysis of the dielectric properties of polar fluids in periodic space, 
Neumann and Steinhauser || derived a measure for the effect of a real-space truncation of Ewald sums, which for 
N = gives Q = [erf(77/2)] 3 . For large 77, the deviation from ideality 5 = 1 — Q 1 ^ 3 scales as 5 ~ exp(— 77 2 /4)/t7, similar 
to what we find based on our analysis of dU/drj. 

The dependence of A on the fc-space cutoff K is depicted in Fig. [l]. The curves obtained from the approximation 



Eq. (12) are compared with the maxima calculated from 10 4 points r randomly chosen in a cubic cell. We observe 
excellent agreement of the approximate formula in the range considered (3 < r/ < 10; < K/2ir < 7). Eq. ( p^ ) 
can therefore be used for assessing the quality of an implementation (N, K, rf) of the Ewald sum. An interesting 
observation from Fig. [l]is that certain K values show a particularly small Fourier-space errors for all 77 values studied. 
The K values are characterized by a gap in the lattice, i.e., there do not exist k vectors such that k 2 = K 2 + 4tt 2 . 
The numerical values are of the form (K/2tt) 2 = 6 + 8m = 6, 14, 22, 30, 38, and 46 0; although followed by a gap, 
(K/2tt) 2 = 27 does not give an optimal cutoff. 

The inverse of Eq. (|9|) can be used to obtain the ratio K/r) given an error in the fc-space sum, 

x{a) P3 2(-lna) 1/2 +ln[4 7 (-ln ( r) 1 / 2 7r- 1 ] (-lna)- 1 / 2 , (13) 

which is asymptotically correct for small errors a, but is already a good approximation for a < 0.1. 

An important observation is that only ratios K/rj enter the formula for the fc-space error Eq. (||) . Correspondingly, 
to achieve the same accuracy with two values 771 and 772, the fc-space cutoff distances have to be chosen proportionally, 
K1/K2 = 771/772- The number of k vectors v(K) for a given value of K scales as K 3 . Thus, the number of k vectors 
required to maintain a given accuracy increases with the third power of the 77 ratios when increasing rj, 

i/(ffi)/i/(jr a ) = {m/m) 3 ■ (u) 

From the analysis of the k-dependent dielectric constant, Neumann || proposed a measure p = 3 exp(— K 2 /4rj 2 ) 
for the Fourier-space error, which also depends only on K/r) and agrees closely with the asymptotic behavior of 
A fe - exp(-X 2 /477 2 )X/7 ? . 

We now determine a maximum useful value of K, given = and 77. This is obtained from a relation K(rj), 
for which the errors of fc-space and real-space truncations are equal, such that a further increase in the number of 
k vectors would not significantly reduce the overall error. Equating the expressions for the real- and Fourier-space 
errors in a cubic lattice, a(K/rj) = 2tt~ 1 I 2 exp(— ?7 2 /4), we obtain an approximate expression, 

K( V ) » ?7 2 + ln(77 2 7 2 /7r) , (15) 

asymptotically valid for large 77. (The relative errors of Eq. (|l|) are less than 0.05 and 0.01 for 77 = 3 and 5, 
respectively.) 

In many calculations, a spherical real-space cutoff R c is introduced in Eq. (|2[), i.e., the argument of the n sum is 
multiplied with a unit step function Q(R C — |r + n|). To find an approximate expression for the numerical error of 
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an implementation (N — 0, 77, K, R c ) where R c is smaller than half of the box length, we use Eq. @ modified by a 
function. We approximate the additional real-space contributions to A as 2ir~ 1 / 2 exp(— i? 2 ?7 2 ), which yields 

A(t7, N = 0, K, R c < 0.5) w (7(^/77) + 2tt~ 1/2 [exp(-?7 2 /4) +exp(-i? 2 ?7 2 )] , (16) 

analogous to Eq. (U). Typically, 77 is large and i? c is chosen smaller than 0.5, such that the exp(— i? 2 7y 2 ) term 
dominates. We can then invert Eq. (|16| ) to find a generalization of Eq. (|l^). This gives the /c-space cutoff K at which 
real- and Fourier-space errors are approximately equal, 

K( V ,R C ) « 2i? c 7 ? 2 + i?- 1 ln(2777i? c 7r- 1 / 2 ) . (17) 

For 77 = 10, we find good agreement for 0.2 < R c < 0.5, with the relative error of Eq. jl^ ) smaller than 2%. 

3 Illustrative examples 

We illustrate our error analysis of Ewald sums using a study of Kusalik Q, who reports relative errors of dipole- 
dipole energies for several configurations of a dipolar soft-sphere fluid calculated with N = 0, 4 < 77 < 7.5, and 
(K/2tt) 2 — 22,30, and 42. Fig. 14 of ref. || shows the relative errors (including the sign) for the three K values as a 
function of 77. Given K, the relative errors are minimal for some values of 77. The optimal combinations of 77 and the 
/c-space cutoff K from Kusalik's calculations are approximately 77 = 5, 5.5, and 6 for (K/2ir) 2 = 22, 30, and 42. These 
values are in excellent agreement with those derived from our analysis, with Eq. ([l5| ) giving 77 = 5.2,5.6, and 6.1 for 
Kusalik's K values. 

In an extension of this study, Kusalik ||Tof reported electrostatic energy, pressure, and dielectric constant of a 
dipolar soft-sphere system for various Ewald-summation parameters (77 = 10, N = 0, 0.236 < R c < 0.31, 46 < 
(K/2ir) 2 < 82). However, the statistical errors-although small-do not allow to establish a conclusive picture, since 
all data are approximately within two estimated standard deviations. Large statistical uncertainties of the order of 
10-20% also prohibit a detailed examination of the errors of the dielectric constant of SPC water, as calculated by 
Belhadj et al. @. 

To further demonstrate the quantitative power of the proposed error analysis, we have studied the energies of 
random configurations of m charges ±1 in a cubic box (with net charge 0). Energies have been calculated for 10 
configurations with m — 8 and m = 32 point charges using Ewald summation (N = 0, K, 77) and by explicit lattice 
sums using lattice vectors n with |n| < 50 (with the correction for a net dipole moment of the box considered ||). 
The relative errors p in the energy with respect to the lattice sums have been determined for (K/2tt) 2 = 10, 20, . . . , 90 
and 3 < 77 < 10. For given values of K, we determine the screening parameter 77 such that the relative errors p(K,r]) 
assume a minimum. 

Fig. ^| shows the relation between optimal K and 77 values together with the derived curve K(rj) fromEq. (|l|). We 
observe excellent qualitative agreement between the derived relation -^(77) and the observed minima. Quantitatively, 
the results for the random configurations suggest somewhat larger /c-space cutoff distances K for given 77. However, 
in most practical applications the charges are more effectively screened by neighboring charges than in random 
configurations, such that the /c-space contributions to the energy tend to be smaller, justifying somewhat smaller K 
values. 

Another important point is the relation between A and the relative errors p(K,rf) in the energy. For p(K,r]), the 
results for the random configurations have been used. A has been calculated from Eqs. ( |l2| ) and (|l5|). Fig. || shows 
minimum values (for given K) of p(K , 77) and A as a function of K. p and A closely follow each other, supporting 
the present error analysis. For the random configurations, they are proportionally related with a factor of about 100. 

4 Discussion 

The present error analysis has important implications on the choice of the Ewald-sum parameters 77, K, and R c in 
computer simulations of condensed-matter systems, helping to avoid unnecessary computational effort and minimize 
the numerical error. The analysis of truncation errors allows to choose 77, K, and R c on a rational basis. Using the 
accuracy measure A, it becomes possible (i) to assess the numerical quality of an Ewald-sum implementation and (ii) 
to compare different implementations using different parameters. 

An important application of the Ewald-summation error analysis in computer-simulation studies is to optimize the 
choice of the screening parameter 77, the real-space cutoff i? c , and the Fourier-space cutoff K regarding computational 



4 



Figure 1: Error A of an Ewald-sum implementation (N = 0, i], K) on a logarithmic scale as a function of the fc-space 
cutoff K/2n. The symbols show the maximum value of A [Eq. (^[)] in the cell estimated from 10 4 random points. The 
lines are calculated from the approximate expression Eq. ([ijj). Results for 3 < t] < 10 are shown. 

speed [m| . Using a few typical configurations of the system, one can minimize the computer time for the 

energy (or force) calculation using combinations of rj and K that give the same error A. The inversion of Eqs. ( |T^ ) 
and (|l6|) yields the appropriate expressions for K , 



where tr _1 is the inverse of cr, as defined in Eq. (||). An approximate analytical expression for a --1 is given in Eq. jl^). 
This strategy of optimizing K) pairs along curves of constant error is particularly important in computational 
studies of large Coulombic systems (e.g., biomolecules in solution), where an overwhelming amount of computer time 
is spent on the calculation of long-range charge interactions. 
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Figure 2: Optimal combination of screening parameter r\ and fc-space cutoff K. The symbols indicate average values 
of r\ that minimize the relative error p(K, rf) in the energy. p{K, rf) has been calculated for 10 random configurations 
of m = 8 (o) and m — 32 (+) charges ±1 in periodically-replicated cubic boxes. The solid line is the predicted result 
from Eq. (15). 



Figure 3: Relative errors p(K , rf) of the energy and Ewald-sum accuracy measure A for optimal combinations of 
screening parameter r\ and fc-space cutoff K. Relative errors p(K , 77) have been calculated for random configurations 
of m = 8 (o) and m = 32 (+) charges ±1; they are compared with A ( ) calculated from Eqs. ( p"2"| ) and flTsj). 
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